Disorder Driven Critical Behavior of Periodic Elastic Media in a Crystal Potential 
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We study a lattice model of a three-dimensional periodic elastic medium at zero temperature with 
exact combinatorial optimization methods. A competition between pinning of the elastic medium, 
representing magnetic flux lines in the mixed phase of a superconductor or charge density waves in a 
crystal, by randomly distributed impurities and a periodic lattice potential gives rise to a continuous 
phase transition from a flat phase to a rough phase. We determine the critical exponents of this 
roughening transition via finite size scaling obtaining v ~ 1.3, /3 ~ 0.05, y/v ~ 2.9 and find that 
they are universal with respect to the periodicity of the lattice potential. The small order parameter 
exponent is reminiscent of the random field Ising critical behavior in 3d. 



A number of materials possess an instability towards 
the formation of a periodically modulated structure in 
space below a particular temperature. Prominent exam- 
ples are charge-density wave systems [[jj , where a Peierls 
instability leads to a state with periodically varying 
charge modulation, or magnetic flux-lines in the mixed 
phase of high-temperature superconductors [0 , where the 
long-range interaction among the lines results in the for- 
mation of the Abrikosov flux-line lattice. Other systems 
forming such periodic structures are spin density waves 
in anisotropic metals, polarization density waves in in- 
commensurate ferroelectrics, and mass density waves in 
superionic conductors [||. 

Usually the periodicity q of this state of broken transla- 
tions! invariance, is incommensurate with the underlying 
crystal lattice, but if q and a reciprocal lattice vector k 
become commensurate (q « k/p) the density wave locks 
in at this wave vector. If fluctuations — either ther- 
mal or induced by impurities, i.e., quenched disorder — 
are weak, these systems are then in a flat phase that 
can be quantified by an order parameter reflecting the 
broken symmetry. Deviations of the local density from 
the perfect periodic structure can be measured by a dis- 
placement field </>(r), which shows long range order in 
the flat phase. When the fluctuations become too strong 
this long-range order vanishes at a roughening transi- 
tion and the system enters a rough phase in which the 
displacement-displacement correlations ([<j>(r) — 0(O)] 2 ) 
diverge with the distance r. 

In the presence of quenched disorder this roughening 
transition is driven by the the competition between two 
pinning forces acting on the periodically modulated flat 
phase: one coming from the underlying lattice potential 
preferring long-range order, the other by the point im- 
purities tending to destroy it. The universality class of 
this transition in the experimentally most relevant case 
of three space dimensions (3d) , its critical exponents and 
the scaling laws, have not been directly scrutinized up 
to now and are the topic of this paper. Thermal fluc- 
tuations are expected to be (dangerously) irrelevant at 
the roughening transition H, and the critical behavior 



at the transition should be dominated by a zero tem- 
perature fixed point analogous to random field critical 
behavior || . Hence the universal properties of the rough- 
ening transition at finite temperatures are expected to be 
identical to the one at zero temperature and the critical 
exponents can in principle be extracted numerically by 
exact ground state calculations |6|,|7J , which is the method 
that we will use here. 

The model Hamiltonian that captures the universal 
properties of the roughening transition under consider- 
ation should contain the following features: It should 
be formulated in terms of a (scalar) displacement field 
4>(r) £ (— oo,+oo); an elastic energy term ^(V</>) 2 as 
the first order (elastic) approximation of the interac- 
tion energy arising from small deformations of the flat 
state <j)(r) = const.; a periodic potential V pcr (<j)) = 
V peT ((j) + 2ir/p) modeling the crystal lattice; and a ran- 
dom potential V ran d{4>) mimicking the effect of impuri- 
ties, which should be invariant under the global shift of 
the whole displacement field <f> — -> <f> + 2-K. The commen- 
surability parameter p entering the periodic potential is 
integer for the lock-in state and is given by the ratio 
of lattice constant of the elastic media with respect to 
that of the underlying periodic potential. The following 
Hamiltonian fulfills these requirements 



H= d d 



7 



| V4>\ 2 — vcos(pcf>) + r/COs(<j) — ip) (1) 



where ip(r) are independent quenched random variables 
uniformly distributed on [— w, 7r] and 7, v, and 7/(r) de- 
note the elastic constant, the periodic potential strength, 
and the random potential strength, respectively. The un- 
derlying elastic approximation for this model is valid as 
long as disorder induced topological defects do not pro- 
liferate. In 2d this actually happens B, but in 3d the 
clastic medium is stable for weak disorder B. 

For v — 0, the Gaussian variational and the functional 
renormalization group (FRG) calculations 1lQ,|lTj] and nu- 
merical studies p2] show that the system is in the elastic 
glass phase, corresponding to a zero-temperature fixed 
point, at all temperatures. The elastic glass phase in 3d 
is characterized by diverging fluctuations 



G(r) = ([0(r o + r)) - 0(r o )] 2 ) ~ 2Aln |r| 



(2) 



at large distances with a universal coefficient A. The 
overbar denotes the disorder average and (. . .) the spatial 
average over ro and the thermal average. 

A simple scaling argument pHl shows that for 
d > 2 the flat phase {4> — 2itn/p, with n a fixed integer) is 
stable as long as the disorder is weak enough: For vanish- 
ing disorder 77 = an excitation (f> — > <fi + 2tt/p over a ter- 
race of linear scale £ costs an elastic energy of the order of 
£ *, whereas for non- vanishing disorder the same exci- 
tation could gain energy of order £ d / 2 . Thus for d > 2 the 
elastic energy loss will dominate over weak disorder and 
the ground state stays flat. Only a strong enough disor- 
der will drive the periodic medium into the rough phase. 
This disorder driven roughening transition was first fig ] 
studied within a variational theory in ||lq| , where a first 
order transition was found, whereas the FRG method 
used in [|4| predicted a continuous roughening transition 
for p > p c (d) = 6/(iry/e) with e = 4 — d at finite disorder 
strength that is determined by a zero-temperature fixed 
point. The order parameter exponent [3 and the correla- 
tion length exponent v were given to leading order in a 
double expansion in e and fj, = p fp\ — 1 by 



4 M 



[3/v= (7r 2 /18)e 



(3) 



A naive insertion of d — 3 and p c = 6/77 into these ex- 
pressions yields values for (3 and v that are incompatible 
with our results which we report now. 

We consider a discrete model for the continuum Hamil- 
tonian (pi). Due to the periodic potential the elastic 
medium remains flat on a microscopic length scale with 



(/) = (2tt/p) h 



h integer . 



(4) 



Therefore, on a coarse-grained level, the medium can be 
described by this integer height variable {h x } represent- 
ing a (3 + l)-dimensional surface on a simple cubic lattice 
with sites x e {1, . . . , L} 3 . Creating steps costs elastic 
energy and the surface is subjected to a random pinning 
potential. These effects plus the periodic potential are 
incorporated in the following solid-on-solid (SOS) model 
Hamiltonian 



^ — 2-^ J(h x ,x);(h y ,y)\h x - h y 



(x,y) 



^ V* cos[(27r/p)/i x - (p x ] 



(5) 



where the first sum runs over nearest neighbor pairs 
(x, y) on a simple cubic lattice. The Hamiltonian has to 
be be invariant under a global shift h — > h + p, which is 
inherited from the symmetry under cf> — > c/)-\-2it of the con- 
tinuum Hamiltonian (pi). Hence we impose a periodicity 
in the step energies J by J(h+ P ,x);(h'+ P ,y) = J(h,*y,(.h',y)- 
Although various aspects of the microscopic physics of 
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FIG. 1. Magnetizations m Ptq (q = 2, 3, 4) at p = 2 and 
A = 2.17,2.20, and 2.23 from top to bottom. The data at 
the critical point A c = 2.20 are connected by solid lines. The 
inset shows the effective exponents for m2,2 (Eq. ([/]))• 

the lattice model (||) and the continuum Hamiltonian (|l|) 
might be different, we can expect the roughening transi- 
tion occurring in both models to be in the same univer- 
sality class, since both models have identical symmetries. 
This 3d SOS model is then mapped onto a (3+l)d ferro- 
magnetic random bond Ising model with an anti-periodic 
boundary condition in the extra-direction, denoted by u 
]l7t . In the ground state the latter induces a 3d inter- 
face, identical to the surface we are searching (up to a 
global shift) that can be determined exactly by using a 
min-cut/max-flow algorithm (see [p|Jl2]| for details). 

We study the model (||) with p — 2,3, and 4 in finite 
L 3 x U lattices with L < 32. The lattice size in the u 
direction, U , has to be chosen large enough in order to 
avoid an interference of the surface with the boundary. 
Random couplings are assigned to the bonds in a unit 
cell of size L 3 x p and they are repeated periodically in 
the u direction. Then the exact ground state configura- 
tion is calculated using a max-flow algorithm. We present 
numerical results obtained by using the uniform distribu- 
tion for < </? x < 2-7T and < i] x < V and the exponential 
distribution, P(J) = J ~ 1 e- J/Jo , for J > 0. The results 
we report do not depend on the choice of the distribution. 
The strength of the random pinning potential is denoted 
by A = V/Jq and we will vary this quantity to trigger 
the roughening transition in our system. 

For each p, we measure the magnetizations 



m M (L,A) = |(e^/ 9 )| (q = 2,3,...) 



(6) 



with ((• • •)) and (• • •) denoting the spatial and the disor- 
der average, respectively, in the ground state. Typically 
the disorder average is taken over 10000 ~ 3000 samples 
for L = 4 ~ 32. Note that the order parameter m = (e 1 ^) 
considered in Ref. B corresponds to m p , q=p , cf., Eq. (H). 
In Fig. pi we show the magnetization for p — 2 as a func- 
tion of L, which scales at the critical point A = A c like 



TABLE I. Estimates for the critical exponents for different 
commensurability parameter p obtained via finite size scaling 
from the numerical data. 



p = 2 
p = 3 
p = 4 



A c 



Ppv/v Pp.a/f PvaI v 



2.20(3) 0.046 (5) 0.034 (3) 

2.475(25) 0.049 (7) 0.037 (9) 

2.95(5) 0.044 (5) 0.033 (5) 



0.022 (3) 1.25(5) 
0.024 (4) 1.29(5) 
0.022 (5) 1.28(8) 



m p,q ~ L~Pv-i' u , with (3 pq and v the order parameter and 
correlation length exponent, respectively. This scaling is 
followed best by the data at A = 2.20, whereas there is a 
downward (upward) curvature for A = 2.23 (2.17) when 
plotting lnm Pi g vs. InL. The critical point A c can be 
directly determined by looking at the effective exponent 



, R i i _ ln(TO p , g (2L)/m p , g (L)) 



(7) 



which is (asymptotically) independent of system size at 
the critical point and equal to the critical exponents 
[3 P ,q/v (see Table I). 

The correlation length exponent is obtained from the 
scaling behavior near the critical point. Each quantity is 
a function of L/£ with the correlation length £ ~ |A — 
Acl - " such that the scaling form of the magnetization is 



•p.q 



(L, A) = L-A>.«/" T((A - Ac)! 1 /") 



(8) 



with a scaling function T . Using the values of A c and 
$i,i I v estimated previously, we determine the correla- 
tion length exponent as the optimal value which yields 
the best data collapse of m2, q =i(L,A). The estimated 
correlation length exponent is also listed in Table I and 
the scaling plot is given in Fig. g (a). 

The correlation length exponent is also determined 
from the susceptibility defined as 



Xp = L 3 (j(e 2 ™'Wp)|2 _ \( e iirih x /p)\ 2 \ 



(9) 



Near the transition point it develops a peak, whose po- 
sition scales as (A*(i) — A c ) ~ L x l v and whose height 
as X*(L) ~ i^l v with the susceptibility exponent 7. For 
each L, A* and x* are obtained by fitting the suscepti- 
bility curve near the peak with a quadratic function, and 
then the critical exponents are extracted to yield that 
i;- 1 = 0.76(5) and 7/1/ = 2.90(5) for p = 2. Both esti- 
mates of v from the magnetization and the susceptibility 
are consistent with each other, and the susceptibility ex- 
ponent satisfies the scaling relation, 7/t/ = d — 1$i,il v 
within the error bars. Fig. @ (b) shows the scaling plot of 
XiL l/u versus (A - A C )L^ V with A c = 2.20, v = 1.25, 
and 7/1/ = 2.90. Except for the smallest system size 
L = 4, those exponents collapse the data well. 

We have performed the same analysis for p = 3 and 
p = 4 and present the critical points and the critical 
exponents in Table I. Fig. (c) and (d) show the corre- 
sponding scaling plots of m p= 3 tq= 3 and m p =4,g=4. 
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FIG. 2. Scaling plots of the magnetization for p = 2 (a), 
p — 3 (c), p = 4 (d), and the susceptibility for p = 2 with 
y/v = 2.90 (b). Parameters used in the plots are listed in 
Table I. 



The order parameter exponents are so small that they 
might also be interpreted as a signature of a first or- 
der transition with /3 = 0. However, they deviate from 
zero systematically for all p and q and, furthermore, the 
scaling analyses show that the transition is associated 
with the diverging correlation length. Therefore we con- 
clude that the transitions are of second order with small 
but finite order parameter exponents, which is reminis- 
cent of the critical behavior of the 3d random-field Ising 
model (RFIM) @,@, for which 0/v = 0.012(4) §. 

Our numerical results deviate from the FRG results |J] 
in many respects. The critical exponents are substan- 
tially different from the analytic results of v ~ 2.59 and 
f3 2 ,i/v ^ 0.548 for p = 2 (see Eq. (g)). Moreover, the or- 
der parameter and the correlation length exponents ap- 
pear to be independent of p within the error bars. 

The discrepancy between FRG and the numerical re- 
sults is surprising since the FRG method produces the 
correct scaling of the elastic glass phase, i.e., the loga- 
rithmic growth of fluctuations as in Eq. (||). Then the 

roughness W 2 



{2K/pf\{hl)-{h^\ ^ExG(x) 
also grows logarithmically in L as W 2 ~ A\nL with 
A = 2jj-(4 — d) |ll|| . These results are confirmed nu- 
merically for a SOS model in 3d @, where W 2 ~ AlnL 
is found with a prefactor A ~ 1.0 in accordance with 
the analytic prediction. We also find this agreement for 
our model (||) in the glass phase (A > A c ): Fig. || (a) 
shows the logarithmic scaling of the roughness W 2 . The 
prefactor A is estimated as \W 2 (2L) - W 2 (L)]/\n2 and 
plotted in Fig. [| (b). The estimates have a strong In- 
dependence. Due to lack of a theory for corrections to 
scaling, we fit them with a quadratic polynomial in 1/ L 
to extrapolate the asymptotic value of A. The results are 
0.98 < A < 1.11, consistent with the analytical and the 
numerical results. If the fluctuations of the elastic glass 
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FIG. 3. (a) W 2 vs. L in semi-log scale, (b) prefactor A of 
W 2 ~ AlnL, (c) effective exponent #3,3 for 7713,3 ~ L~ 03 - 3 , 
and (d) exponent ratios r = #3,2 /03, 3 (empty symbols) and 
#3,2 /#3, 4 (filled symbols). All data are made from p = 3 in 
the disordered phase. 



phase have the Gaussian nature as assumed in Ref. [[Ll) , 
the roughness, W 2 ~ ^41nL, and the magnetizations, 
m p,q ~ L~ 0p - q , are not independent quantities since they 
should obey ( e ^W-^(r')) ^ e -±<W>(r)-</>(r')] 2 > implying 
the relation 



A = 2 6. 



p.p 



/Op 



J p,ql w p,q 



q' 2 /q 2 



(10) 



Fig. H (c) shows the effective exponents $3,3. The polyno- 
mial fitting is used to extrapolate the asymptotic value, 
which yields 0.51 < #3,3 < 0.62. We also calculate 
1117713,2/1117713,3 and In 7773, 2/ln 7713,4, which approach 2.25 
and 4.0, respectively, as L increases, cf., Fig. |(d). Those 
values satisfy the scaling relations in Eq. (lCJ) approxi- 
mately. This is a strong evidence for the Gaussian nature 
of the fluctuations in the elastic glass phase and hence 
justifies the analytic approaches in this regime. 

However, when using the same procedure as in the 
glass phase to determine A c at the critical point for 
p = 2,3,4, we obtain A c ~ 0.18 (p = 2), 0.092 (p = 3), 
and 0.046 (p = 4). Note that A c appears to be approx- 
imately inversely proportional to p 2 , which implies that 
the bare width Wq = (/i 2 .) — (h x ) 2 is independent of p. 
The Gaussian theory requires that 



A r 



2(/W*0 . Pp., q /Pp., q ' = q'Vq 2 



(11) 



but A c and [fip V jv) listed in Table I violate the first 
relation by a factor of two. The ratios P p ^/Pp,3 — 1-5 
and fip^/ PpA — 2-3 are also far from the values 9/4 and 
4, required by ([ll]). This implies a strongly non-Gaussian 
nature of the fixed point in 3d and provides a hint why the 
FRG-prediction (pf) for 3d differs from ours. On the other 
hand, the latter are based on a double expansion around 
d = 4 and p = p c and it is well possible that d = 3 and the 
values for p we have considered here are simply beyond 



the validity of such first order perturbation expansion. 
We want to stress that we think that the system sizes we 
studied are sufficiently large to see the true asymptotic 
behavior of the roughening transition since we capture 
correctly the features of the rough phase fixed point. 

In summary we presented the first numerical study of a 
disorder driven roughening transition in a periodic elastic 
medium. Our results for the critical exponents deviate 
significantly from the predictions of a recent analytical 
FRG calculation and we discussed the validity of its un- 
derlying Gaussian approximation. We found that this 
new universality class is reminiscent of random field crit- 
ical behavior in 3d including a very small order parameter 
exponent. To complete the picture of the underlying zero 
temperature fixed point scenario, one has to compute the 
violation of hyperscaling exponent 9, which necessitates 
techniques different from those used in this work p9[. 
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